package com.appspot.seqpop.service;

import java.io.Serializable;

public class Coefficient implements Serializable {

    private static final long serialVersionUID = 1L;

    private int sampleCount;
    
    private double a;
    private double b;
    private double e1 = Double.NaN;
    private double e2 = Double.NaN;
    private double udStar = Double.NaN;
    private double vdStar = Double.NaN;
    private double ufStar = Double.NaN;
    private double vfStar = Double.NaN;
    
    public Coefficient(int n) {
        if (n < 2) {
            throw new IllegalArgumentException();
        }
        this.sampleCount = n;
        a = computeA(n);
        b = computeB(n);
//        double d = computeD(n, a);
//        e1 = computeE1(n, a);
//        e2 = computeE2(n, a, b);
//        vdStar = computeVDStar(n, a, b, d);
//        udStar = computeUDStar(n, a, vdStar);
//        vfStar = computeVFStar(n, a, b);
//        ufStar = computeUFStar(n, a, vfStar);
    }
    
    private static double computeA(final int n) {
        double sum = 0.0;
        for (int i = 1; i < n; i++) {
            sum += 1.0 / i;
        }
        return sum;
    }
    
    private static double computeB(final int n) {
        double sum = 0.0;
        for (int i = 1; i < n; i++) {
            sum += 1.0 / (i * i);
        }
        return sum;
    }
    
    private static double computeE1(final int n, final double a1) {
        double b1 = (double)(n + 1) / (3 * (n - 1));
        double c1 = b1 - 1 / a1;
        return c1 / a1;
    }
    
    private static double computeE2(final int n, final double a1, final double a2) {
        double b2 = (double)(2 * (n * n + n + 3)) / (9 * n * (n - 1));
        double c2 = b2 - (n + 2) / (a1 * n) + a2 / (a1 * a1);
        return c2 / (a1 * a1 + a2);
    }
    
    private static double computeD(final int n, final double a) {
        double an_1 = a + 1.0 / n;
        double c = 2 * (n * a - 2 * (n - 1)) / ((n - 1) * (n - 2));
        return c + (double)(n - 2) / ((n - 1) * (n - 1)) + 2.0 / (n - 1) * (1.5 - (2 * an_1 - 3) / (n - 2) - 1.0 / n);
    }
    
    private static double computeVDStar(final int n, final double a, final double b, final double d) {
        return ((double)(n * n) / ((n - 1) * (n - 1)) * b + a * a * d - 2 * n * a * (a + 1) / ((n - 1) * (n - 1))) / (a * a + b);
    }
    
    private static double computeUDStar(final int n, final double a, final double v) {
        return (double)n / (n - 1) * (a - (double)n / (n - 1)) - v;
    }
    
    private static double computeVFStar(final int n, final double a, final double b) {
        double v1 = (double)(2 * n * n * n + 110 * n * n - 255 * n + 153) / (9 * n * n * (n - 1));
        double v2 = 2 * (n - 1) * a / (n * n);
        double v3 = 8 * b / n;
        return (v1 + v2 - v3) / (a * a + b);
    }
    
    private static double computeUFStar(final int n, final double a, final double v) {
        double an_1 = a + 1.0 / n;
        return (4 * n * n + 19 * n + 3 - 12 * (n + 1) * an_1) / (3 * n * (n - 1)) / a - v;
    }

    public int getSampleCount() {
        return this.sampleCount;
    }

    /**
     * The coefficient of Watterson's theta estimator (Watterson 1975).
     * Let a and S be the coefficient and the number of segregating sites,
     * then E[θ] = a*S.
     * @param n The number of sample sequences.
     * @return The a described above.
     */
    public double getA() {
        return this.a;
    }

    /**
     * The coefficient of Watterson's theta estimator (Watterson 1975).
     * Let a and b be the coefficients from {@linkplain #getA()} and {@linkplain #getB()},
     * and let S be the number of segregating sites,
     * then Var[θ] = a*S + b*S^2.
     * @param n The number of sample sequences.
     * @return The b described above.
     */
    public double getB() {
        return this.b;
    }

    /**
     * The coefficient of the variance of Tajima's D statistic (Tajima 1989).
     * Let S and Π be the number of segregating sites and the average number of
     * nucleotide differences, and let a, e1 and e2 be the coefficients from
     * {@linkplain #getA()}, {@linkplain #getE1()} and {@linkplain #getE2()},
     * then Var[Π - S/a] = e1*S + e2*S(S-1) 
     * @return The e1 described above.
     */
    public double getE1() {
        if (Double.isNaN(this.e1)) {
            synchronized (this) {
                if (Double.isNaN(this.e1)) {
                    int n = getSampleCount();
                    double a = getA();
                    this.e1 = computeE1(n, a);
                }
            }
        }
        return this.e1;
    }

    /**
     * The coefficient of the variance of Tajima's D statistic (Tajima 1989).
     * Let S and Π be the number of segregating sites and the average number of
     * nucleotide differences, and let a, e1 and e2 be the coefficients from
     * {@linkplain #getA()}, {@linkplain #getE1()} and {@linkplain #getE2()},
     * then Var[Π - S/a] = e1*S + e2*S(S-1) 
     * @return The e2 described above.
     */
    public double getE2() {
        if (Double.isNaN(this.e2)) {
            synchronized (this) {
                if (Double.isNaN(this.e2)) {
                    int n = getSampleCount();
                    double a = getA();
                    double b = getB();
                    this.e2 = computeE2(n, a, b);
                }
            }
        }
        return this.e2;
    }

    public double getUDStar() {
        if (Double.isNaN(this.udStar)) {
            synchronized (this) {
                if (Double.isNaN(this.udStar)) {
                    int n = getSampleCount();
                    double a = getA();
                    double v = getVDStar();
                    this.udStar = computeUDStar(n, a, v);
                }
            }
        }
        return this.udStar;
    }

    public double getVDStar() {
        if (Double.isNaN(this.vdStar)) {
            synchronized (this) {
                if (Double.isNaN(this.vdStar)) {
                    int n = getSampleCount();
                    double a = getA();
                    double b = getB();
                    double d = computeD(n, a);
                    this.vdStar = computeVDStar(n, a, b, d);
                }
            }
        }
        return this.vdStar;
    }

    public double getUFStar() {
        if (Double.isNaN(this.ufStar)) {
            synchronized (this) {
                if (Double.isNaN(this.ufStar)) {
                    int n = getSampleCount();
                    double a = getA();
                    double v = getVFStar();
                    this.ufStar = computeUFStar(n, a, v);
                }
            }
        }
        return this.ufStar;
    }

    public double getVFStar() {
        if (Double.isNaN(this.vfStar)) {
            synchronized (this) {
                if (Double.isNaN(this.vfStar)) {
                    int n = getSampleCount();
                    double a = getA();
                    double b = getB();
                    this.vfStar = computeVFStar(n, a, b);
                }
            }
        }
        return this.vfStar;
    }
}
